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ABSTRACT 

We study the influence of the choice of cooling function on the formation of molec- 
ular clouds in high-resolution three-dimensional simulations of converging flows. We 
directly comp are the results obtained usin g the simple, parametrized cooling function 
introduced by iKovama Sz Inutsukal (|2002l) and used by a number of converging flow 
studies with the results of the detailed calculation of the non-equilibrium chemistry 
and thermal balance of the gas. We find that a number of the cloud properties, such 
as the mass and volume filling fractions of cold gas, are relatively insensitive to the 
choice of cooling function. On the other hand, the cloud morphology and the large- 
scale velocity distribution of the gas do strongly depend on the cooling function. We 
show that the differences that we see can largely be explained by differences in the 
way that Lyman-a cooling is treated in the two complementary approaches, and that 
a proper non-equilibrium treatment of the ionisation and recombination of the gas is 
necessary in order to model the high-temperature cooling correctly. 

We also investigate the properties of the dense clumps formed within the cloud. 
In agreement with previous models, we find that the majority of these clumps are not 
self-gravitating, suggesting that some form of large-scale collapse of the cloud may be 
required in order to produce gravitationally unstable clumps and hence stars. Overall, 
the physical properties of the dense clumps are similar in both simulations, suggesting 
that they do not depend strongly on the choice of cooling function. However, we do 
find a systematic difference of around 10 K in the mean temperatures of the clumps 
produced by the two models. 

Key words: astrochemistry - molecular processes - ISM: clouds - ISM: molecules - 
methods: numerical - turbulence. 



1 INTRODUCTION 

An important goal of star formation research is to under- 
stand the formation and evolution of giant molecular clouds 
(GMCs), known to be the hosts of all observed Galactic 
star formation. However, as yet this remains an unsolved 
problem. Traditionally, molecular clouds were viewed as be- 
ing virialized structures in the interstellar medium (ISM), 
having relatively long lifetimes and a significant delay be- 
tween the formati on of the cloud and the onset of star 
formation (see e.g. Zuckerman fc Palmer! Il974l : IWoodwardl 
Il978l : IScoville fc Hershlll979l l. However, recent studies have 
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suggested that GMCs begin forming stars shortly after 
they th emselves form, and are non-e q uilibrium entities 
(see e.g. iBallesteros-Paredes et al.lll999dlbllElmegreen]|200i 
lHartmann et al.ll200ll : Vazquez-Semadeni et al .Il2003l . l200 



Their formation and evolution are dominat ed by the ef- 



fects o f supe r sonic turbulent motions (see e.g. iLarsonl 



Myers! Il983l: ISolomon et al.l 1 19871 : iFalgarone etal. 
Hever fc Brunti2004l and the revie ws bvlMac Low fc K 



20041: lElmegreen fc Scald 12004 IScalo fc Elmegreenl 



1981 



1992 



2004; 



McKee fc Ostrikerl 120071 ) and are rapid, with a timescale 



comparable to those of important chemical processes such 
as the conversion of atomic to molecular hydrogen or the 
freczc-out of molecules onto the surfaces of interstellar dust 
grains. Therefore, the dynamics and chemistry of the gas are 
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strongly coupled, with one directly influencing the evolution 
of the other, and they must be modelled together. 

A promising theory for producing GMCs, as well 
as for generating turbulence within them, suggests that 
these clouds form in places where streams of warm 
atomic gas collide. Work by a number of different 
groups has shown that the dense sheets and fila- 
ments that build up at the interf ace of the colliding 
flows become thermally unstable (e.g. Audit fc H ennebclle 
2005l;IVazquez-Semadeni et ah 2007; Hennebelle et al.ll200sf 



Baneriee et~aH 2009: Hcitsch & Hartmann 2008). The com- 



pressed unstable gas rapidly radiates away most of its ther- 
mal energy, significantly decreasing its temperature. As the 
temperature of the gas falls, it is compressed by the warmer 
material surrounding it, and so the large drop in tempera- 
ture is associated with a large increase in density. In some 
regions, strong pressure gradients are created that then act 
to drive turbulent flows throughout the interaction region. 
The outcome is a set of dense, cold clumps of gas, embed- 
ded in a more diffuse, turbulent flow, and with a filamentary 
morphology reminiscent of that found by observations of real 
GMC s (e.g. iMen'shchikov et al.l |2010| ; lArzoumanian et alj 
2011). Since the dense clumps observed in real GMCs are 
known to be the sites at which stars form, these results sug- 
gest that there is an important link between the large scale 
physics of GMC formation and the smaller scale physics of 
star formation. 

Unfortunately, most three-dimensional studies of this 
process performed so far make use of highly simplified treat- 
ments of the thermal energy balance of the gas. Heating and 
cooling are modelled using simplified parameterizations that 
assume that the rates are functions only of density and tem- 
perature, and that ignore the effects of chemical changes 
or dust shielding (s ee e.g. IVazquez-Semadeni et alj 120071; 
iBaneriee et al.ll2009h . No effort has been made to determine 
whether this simplified approach is adequate for describing 
the dynamics of the gas, or whether we would expect to find 
significant changes in behaviour if we were to use a more 
detailed chemistry and cooling model. 

In this paper, we attempt to improve on this situation. 
We perform high-resolution 3D simulations of cloud forma- 
tion in colliding flows, and directly compare the results ob- 
tained from a simplified cooling model with those that we 
obtain from a self-consistent treatment of the cooling and 
chemistry of the gas. Our main goal is to understand how 
the use of a more accurate thermal model affects the dy- 
namics of the flow and the nature of the structures that 
form within it. 



2 NUMERICAL MODEL 

2.1 The numerical code and setup 

We consider the atomic phase of the interstellar medium 
(ISM), whose behaviour is governed by the equations 



dp 
dt 



+ V ■ (pv) = 



(1) 
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where p is the gas density, v is the fluid velocity, B is the 
magnetic field strength, which also must satisfy the con- 
straint V B = 0, E = P/(7-l) + pH 2 /2 is the total energy 
per unit volume, P is the thermal pressure, and 7 = 5/3 is 
the adiabatic index. In the energy equation, n is the num- 
ber density of hydrogen nuclei, which is related to the mass 
density via n = p/(1.4m p ), where m p is the proton mass, 
nF is the radiative heating rate per unit volume and n 2 A is 
the radiative cooling rate per unit volume. In writing down 
this set of equations, we have assumed that we can neglect 
the effects of the self-gravity of the gas. We make a similar 
assumption in the simulations presented in this paper, and 
defer investigation of the more computationally demanding 
self-gravitating case to future work. 

We model the collision of two large cylindrical streams 
of warm atomic gas using a modified ve rsion of the adaptiv e 
mesh refinement (AMR) code FLASH (|Frvxell et al.ll200Gh . 
Our modifications include the addition of a simplified but 
accurate treatment of the most important hydrogen chem- 
istry, together with a detailed atomic an d molecular coolin g 
function. They are described in detail in iMicic et al.l l|2012l ). 
For our study we use the MHD 5-wave Bouchut solver 
whic h preserves positive states for the density and pres- 
sure llBouchut. Klingenberg fc Waa gan 2007 1. |2010| ; IWaaganl 
120091 ; IWaagan. Federrath fc Klingenberd l201ll) . This im- 
plementation also applies the truncation-error method 
|Powell et all Fl 999) to maintain small V • B errors and to 
avoid unphysical magnetic tension terms. 

W e use a similar setup to that studied in lBaneriee et all 
(2009), which itself was based on model L256Au0.17 from 
the study of lVazquez-Semadeni et all (|2007l ) . The two cylin- 
drical streams, each 112 pc long and 32 pc in radius, are 
given an initial, slightly supersonic inflow velocity so that 
they collide at the centre of the numerical box (x = pc). 
The (256 pc) 3 simulation box is periodic, and the streams 
are completely contained within it, such that the resulting 
cloud occupies a relatively small volume far from the bound- 
aries, and interacts freely with its diffuse environment, with 
relatively little effect from the boundaries. The box is ini- 
tially filled with warm atomic gas with a uniform number 
density n = 1 cm" 3 . This corresponds to an initial mass 
density of p = 2.12xl0~ 24 g cm~ 3 , if we adopt a 10:1 ratio 
of hydrogen to helium, by number. The initial temperature 
of the atomic gas is ~ 5000K, corresponding to an isother- 
mal sound speed of 5.7 km s . The initial velocity of each 
flow is 7 km s -1 , and so the flows have an initial isother- 
mal Mach number of 1.22. At temperature of T = 5000K 
for the warm phase, this implies that the cold phase comes 
into hydrostatic thermal pressure balance w ith the warm 
gas a t a density of roughl y 100 cm -3 (see e.g. IWolfire et al.l 
I1995I : IWolfire et ai1l2003T ). Furthermore, we add 10% ran- 
dom velocity perturbations to the bulk stream. Finally, we 
note that we include a magnetic field, which we assume 
to be oriented parallel to the inflow. The initial magnetic 
field strength is taken to be 3 pG, consistent with esti- 
mates of the mean Galactic magnetic field strength (|Beckl 
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200 ll). This corresponds to a c ritical mass-to-flux ratio (see 



Vazguez-Semadeni et aljfeoill ) 



Table 1. Reactions in our non-equilibrium chemical model. 



We follow the collision of the streams of gas with up 
to 11 AMR refinement levels corresponding to a maximum 
effective resolution of 8192 3 grid cells, or a grid spacing of 
Ax = 0.03 pc in each direction. We use a Jeans- type cri- 
terion (|Truelove et al. 1 1 19971 ; iFederrath et~afll201ll ) for the 
dynamical mesh refinement which requires that the local 
Jeans length is resolved with at least 10 grid cells while re- 
finement is active. Our use of a Jeans-type criterion ensures 
that we resolve any structures large enough to be gravita- 
tionally unstable (although as we do not include the effects 
of self-gravity in these simulations, these structures will not 
actually undergo gravitational collapse) . In practice, we find 
that this criterion also allows us to resolve a considerable 
amount of the smaller-scale structure of the g discussed 
in more detail in Section 3. Nevertheless, our limited reso- 
lution means that we will still miss structures on very small 
scales, as we expect the thermal instability to create struc- 
tures on all scales larger than the Field length i|Fieldll 19651 ). 
which in th e CNM is typically of the order of 10~ 3 pc or 
less (see e.g. iGlover fc Mac Lo w 2007a). Fortunately, previ- 
ous studies have shown that the failure to resolve this very 
small scale structure does not significantly affect the dy - 
namics of the gas on larger scales (see e.g. iGressell [20091) . 
and hence this limitation should not significantly affect our 
results. 



2.2 Cooling and heating 

The simulatio ns of [ vkzouez-Semadeni et al.l (|2007h and 
iBaneriee et al.l (|2009h used a cooling function derived 
from the one-dimension al colliding flow models of 
iKovama fc Inut suka ( 200 ^). A simpl e analytical fit was pro- 
vided bv lKovama fc Inutsukal (|2002j ). and has the forrro 



r = 2.0 x 10~ 2t> erg s , 
A(T) - ( 1.184 x IQ 5 \ 

— = exp (r r + iooo ) 

+ 1.4 x 10~ 2 VTexp f-p) cm3 - 



(5) 



(6) 



where T is the gas temperature in Kelvin. With this cool- 
ing function, the simulated ISM is thermally unstable in the 
density range 1 < n < 10 cm -3 , corresponding to equilib- 
rium temperatures in the range 500 < T < 5000 K. 

This simplified treatment of the heating and cooling 
of the gas does not account for any chemical effects. It 
assumes that the photoelectric heating efficiency is con- 
stant, whereas in practice it is kn own to depend on the 
electron number d ensity (see e.g. iBakes fc Tielensl 1 19941 : 
IWolfire et aDll995h . It also assumes that Lyman- a cool- 
ing is extremely efficient: the high tempera ture cooling rate 
assumed in the IKovama fc Inutsukal ([2002J) function corre- 
sponds to the cooling one would expect in a gas in which the 



1 Not e that the version of this fit printed in IKovama fc Inutsukal 
suffers from a significant typographical error, which has the 
effect of making the low-temperature cooling rate far too large. 
The ver sion we quote here is the correct ed version of their fit, as 
given in IVazqucz-Scmadcni ct al. (20q3). 
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Notes: "c.r." denotes a cosmic ray particle, and 7x denotes an 

X-ray photon 
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electron and atomic hydrogen number densities are approx- 
imately equal. In our present study, we investigate how the 
results that we obtain when we properly account for these 
effects compare with the results produced by this more sim- 
plified treatment. 

To do this, we perform two simulations using 



the Kovama & 


nutsukal 


IBaneriee et al.l ( 


20091). I 



dependent chemical model and cooling function described in 
IGlover fc Mac Low! ( 2007a, b). This model follows the abun- 
dances of four chemical species - free electrons, H + , H, and 
H 2 - linked by the reactions listed in Table [1] It assumes 
that any carbon in the gas remains in the form of C + and 
that any oxygen remains in atomic form, and so underesti- 
mate s the cooling rate in re gions dominated by CO. How- 
ever. Gl over fc Clarkl (|2012al ) have shown that this does not 
have a strong effect on the dynamics of the gas on scales 
larger than those of individual pre-stellar cores, and so mak- 
ing this simplifying assumption should not greatly affect our 
results. 

We use an implementation of the Glover & Mac Low 
model within FL ASH that is described in detail in 
iMicic et al l l|2012l ). Our approach uses FLASH'S standard 
tracer field implementation to directly follow the advection 
of the fractional abundances of molecular hydrogen (ih 2 ) 
and ionised hydrogen (x H +). The abundances of the other 
two species - atomic hydrogen (ih) and electrons (x ) - are 
computed from the conservation laws for charge 

x c = x n + + x c + + x Si + (7) 

and for the total amount of hydrogen 

XH = XH,tot — X H + — 2XH 2 ■ (8) 

Here, Krr.tot is the total abundance of hydrogen nuclei in all 
forms, which we normalize to unity, and x c + and x Si + are 
the abundances of ionised carbon and silicon, respectively, 
which remain fixed throughout the simulations. 

We assume that carbon, oxygen and silicon remain in 
the form of C + , O and Si + throughout the simulation, 
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Figure 1. Mass- weighted density (left panel) and temperature (right panel) probability distribution function (PDF) at time t = 22 Myr. 
The solid line presents the PDF in the run with non-equilibrium chemical model, while the dashed line corresponds to the PDF in the 
run with Koyama & Inutsuka cooling function. 



and adopt fractional abundances for these species given by 
x c+ = 1.41 x 10~ 4 , xp = 3.16 x IP " 4 and i Si+ = 1.5 x 1(T 5 , 
respectively l|Sembach et alJfeoOCh . The radiative and chem- 
ical heating and cooling of the gas is modelled with a 
cooling function that contains contributions from a range 
of processes, of which the most important are C + and 
O fine structure cooling, Lyman-a cooling, and photoelec- 
tric heating. Full details of these processes, along with the 
other contributions to our cooling function, can be fo und in 
iGlover fc Mac Low! (|2007al lbf) and lMicic et alj \2Q\t \. 

In our run with the Glover & Mac Low chemistry and 
cooling model, we adopt a value of Go = 1.0 for the strength 
of the interste llar radiation field in units of the Habing field 
ilHabindlilxSsh . We take the cosmic ray ionisation rate of 
atomic hydrogen to be £h = 10 -17 s _1 and assume that the 
ratio of this rate to the cosmic ray ionisation rate of H2 is 
the same as given in the UMIST a strochemistry database 
|Le Teuff. Millar, fc Markwickll20odh . We include the effects 
of X-ray ionisation and heating using t he prescription given 
in Appendix A of Wol fire et al.1 |l995), and assume a uni- 
form absorbing column density of warm atomic hydrogen 
iV w = 10 19 cm -3 . In our present study, we do not include 
the effects of self-gravity, dust shielding, or H2 self-shielding. 
We note that although we expect dust shielding to have a 
significant effect on the thermal state of the gas regions with 
mean visual extinctions Av > 1, these account for only a 
small fraction of the simulation volume. 



3 RESULTS 

3.1 Density and temperature distributions 

The sequence of events that occurs within our two sim- 
ulations is broadly similar in both cases, and is also in 
goo d agreement with the results o f previous studies (se e 
e.g. IVazquez-Semadeni et al' I 120071 ; iBaneriee et al. 2009). 
We therefore begin by briefly describing this sequence of 
events, before moving on to look at the differences that do 
occur between the two runs. 

At the interface where our transonic, converging flows 
collide, the gas is shocked and moderately compressed. This 



compression destabilises the gas, triggering a thermal insta- 
bility (TI) that causes the gas to cool rapidly. As the gas 
cools, it is compressed by the thermal pressure of the sur- 
rounding warm gas, leading to a rapid increase in its den- 
sity. This process comes to an end once the gas reaches the 
equilibrium temperature of the cold neutral medium (CNM) 
phase, which for the conditions simulated here is below 
100 K. The cool dense gas initially forms a sheet that then 
fragments into filaments and ultimately into small, pressure- 
confined clumps. As the thermal pressure of the dense gas is 
in close balance with the total (thermal plus ram) pressure 
of the warm neutral medium (WNM) outside it, the cold gas 
can easily reach number densities of the order of 100 cm -3 , 
compar able to the mean density of the gas in many G MCs 
(sec e.g. Roman- Duval et aUbo'ld : iHughes et aLll2010l ). 

The cloud of gas that forms in the interface region is 
composed of a mixture of diffuse and dense gas, includ- 
ing a significant fraction of material in the thermally un- 
stable region intermediate between the CNM and WNM 
phases (see Figure [T]). Rather than the classical picture of 
a two-phase medium, we find instead a continuous distri- 
bution of densities and temperatures, albeit one with clear 
pea ks corresponding to the CNM and WNM regimes (see 
alsolVazquez-Semadeni et al.ll2000l ; iGazol et al.ll200ll . |2005l : 
I Audit fc Hennebellell2005l ). 

In Figure [T] we plot mass- weighted density (left-hand 
panel) and temperature (right-hand panel) probability dis- 
tribution functions (PDFs) for both models at a time t = 
22 Myr, several million years after the end of the inflow. It 
is clear from the Figure that these PDFs do not differ by 
much between the two runs. The main difference that is ap- 
parent is that the temperatures that we recover for the CNM 
and WNM phases (the two clear peaks in the temperature 
PDF) are slightly smaller in the simulation run using the 
Koyama & Inutsuka cooling function than in the simulation 
using the Glover & Mac Low treatment. This difference in 
behaviour is relatively simple to understand. At high tem- 
peratures (T > 7000 K), the Koyama & Inutsuka cooling 
rate coefficient is given approximately by 



Aki(T) ~ 2 x 10" 



exp 



1.184 x 10 5 
T+1000 



(9) 
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Figure 2. Two-dimensional PDFs of temperature and density for the two simulations, for a time t = 22 Myr. The left-hand panel shows 
the results from the simulation that used the Koyama & Inutsuka cooling function, and the right-hand panel shows the results from the 
simulation using the full non-equilibrium treatment. The fraction of mass in each region of the density- temperature space is indicated 
by the colour scale. The solid line shows the equilibrium temperature as a function of density, derived under the assumption that the 
gas is also in chemical equilibrium. 



The main coolant in the Glover & Mac Low treatment at 
these temperatures is Lyman -Q cooling, for which they use 
the following expression from lCenl l|l992h 



A. 



7.5 x 10" 



Ly-a = — ; , = exp — XeXH, (10) 

where x a = n e /n is the fractional abundance of electrons 
and xh = rm/n is the fractional abundance of atomic hy- 
drogen. Comparing these two cooling rates, we find that 
they produce comparable amounts of cooling only when 
x c ~ ih — 1/2, i.e. only when the chemical state of the gas is 
such that we get roughly the maximum amount of Lyman- 
a cooling possible. In highly ionised gas with n a 2> uh, 
or predominantly neutral gas with n c <C nu, the Lyman- 
a cooling rate is considerably smaller, and hence in these 
conditions, the Koyama & Inutsuka treatment significantly 
overestimates the true cooling rate of the gas. The expo- 
nential temperature dependence of the cooling rate means 
that a large error in the value of the rate leads to only a 
small error in the gas temperature, but this is sufficient to 
explain the offset in the characteristic temperature of the 
WNM that we find when we compare our two simulations. 

The difference in the CNM temperatures is not caused 
by a difference in the low temperature cooling rates, but 
rather by a difference in the radiative heating rate. In the 
Koyama & Inutsuka treatment, the heating rate throughout 
the gas is simply T = 2 x 10~ 26 erg s _1 . In the Glover & 
Mac Low treatment, on the other hand, the heating rate 
is sensitive to the chemical composition of the gas ow- 
ing to the fact that the photoelect ric heating efficiency is 
a function of the electron de nsity (|Bakes fc Tielena 1 19941 ; 
IWeingartner fc Draind l2001bh . In cold, dense gas, photo- 
electric heating is relatively efficient, and the heating rate 
is given approximately by r pc ~ 5 x 10 -26 erg s _1 , i.e. it 
is roughly 2.5 times larger than assumed in the Koyama & 
Inutsuka treatment. It is therefore not surprising that the 
gas can cool to somewhat lower temperatures in this case. 



In order to investigate the dependence of temperature 
on density in our two simulations, we plot in Figure [2] the 
two-dimensional (2D) PDFs of these quantities. The results 
for the Koyama & Inutsuka run are shown in the left-hand 
panel, while the results of the non-equilibrium run are shown 
in the right-hand panel. The equilibrium temperature of the 
gas is indicated using the solid line. In the case of the sim- 
ulation using the non-equilibrium treatment of chemistry 
and cooling, we derived an equilibrium temperature at each 
density by assuming that the gas was also in chemical equi- 
librium. We see from Figure [5] that in the simulation with 
the Koyama & Inutsuka cooling function, most of the gas 
has a temperature close to the equilibrium value. In the 
simulation with the Glover & Mac Low chemistry and cool- 
ing, on the other hand, the departures from equilibrium are 
more pronounced. The gas is close to the equilibrium tem- 
perature at densities n < 1 cm -3 and n > 30 cm -3 , but the 
equilibrium temperature curve does not give a good descrip- 
tion of the distribution of gas temperatures at intermediate 
densities. In this intermediate regime, the temperature falls 
off less rapidly with increasing density than predicted by 
the equilibrium temperature curve. This is a consequence 
of the sensitivity of the photoelectric heating rate to the 
electron number density, n c . As n e increases, the net posi- 
tive charge of the dust grain population decreases, or even 
becomes negative. This makes it easier for incoming pho- 
tons to cause the ejection of photoelectrons. Consequently, 
the photoelectric heating efficiency tends to increase with 
increasing electr on number density, up to a limiting value 
of a few percent (iBakes fc Tielendll993 ; IWolfire et aLlll995l ; 



IWeingartner fc Draineil2001bl ). In the density and tempera- 
ture regime where we see the greatest deviations from the 
equilibrium temperature curve, the photoelectric heating ef- 
ficiency has not yet reached this limiting value, and hence 
any difference between the actual electron number density 
and the equilibrium value leads to a photoelectric heating 
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Figure 3. Evolution with time in our two simulations of the mass (left-hand panel) and volume (right-hand panel) fractions of cold gas, 
defined here as gas with temperature T ^ 300K. The solid line is for the run with the non-equilibrium chemical model, while the dashed 
line is for the run with the Koyama & Inutsuka cooling function. 



rate that also differs from the value that it would have in 
equilibrium. As the recombination timescale in this gas is 
relatively long, of the order of a few Myr or more, the gas is 
generally slightly more ionised than it would be in equilib- 
rium, and hence is heated slightly more efficiently. 

Finally, as we know that stars form in cold gas, it is 
interesting to examine how the cold gas fraction differs in our 
two simulations. As we have already seen, the cold and warm 
phases in our simulations are not completely distinct, and a 
significant fraction of the mass of the cloud lies intermediate 
between these two phases. The definition of "cold" gas is 
therefore somewhat subjective. In this paper, we define cold 
gas to be gas with a temperature T < 300 K. (We note 
from Figure [2] that in both simulations, the majority of the 
gas with a temperature as low as this has a density n > 
10 cm -3 ). In Figure [3] we plot how the mass fraction (left- 
hand panel) and volume fraction (right-hand panel) of the 
cold gas evolves with time in both simulations. 

Prior to the end of the inflow, at t ~ 16 Myr, both sim- 
ulations show extremely similar behaviour. The amount of 
cold gas is small - it accounts for only 5% of the total mass 
and only a tiny fraction of the total volume of the simulation 
- and the cold gas fraction increases only slowly with time. 
After the end of the inflow, however, greater differences be- 
come apparent between the two runs. In the run with the 
non-equilibrium cooling and chemistry, the cold gas fraction 
remains small up to t ~ 25 Myr, but thereafter begins to 
increase rapidly. On the other hand, in the run with the 
Koyama & Inutsuka cooling function, the cold gas fraction 
grows steadily from time t ~ 20 Myr until the end of the 
simulation, at a somewhat faster rate than during the in- 
flow phase. In the interval 20 < t < 28 Myr, there is more 
cold gas in the Koyama & Inutsuka run than in the non- 
equilibrium run, but at later times the latter run has the 
most cold gas. Figure [3] also indicates that the volume fill- 
ing factor of the cold gas in the Koyama & Inutsuka run is 
larger than in the other run, although the difference between 
the two is not large. We can understand the difference in the 
behaviour of the cold gas fraction in the two run by look- 
ing at how the velocity field generated in the cloud differs 
between the runs, which we examine in the next section. 



3.2 Cloud structure and velocity 

In Figure [4] we show projections of the column density of 
hydrogen nuclei, N, along the axis of the flow at three com- 
parable output times for both of our simulations. The top 
panels show the results for the run with the Koyama & In- 
utsuka cooling function, while the bottom panels show the 
results from our non-equilibrium chemistry run. Figure [5] 
shows a similar comparison, but for a direction perpendicu- 
lar to the flow. 

The images in Figure U show us clearly that the cloud 
is not a homogeneous entity, but rather is composed of nu- 
merous dense clumps embedded in lower density filaments. 
Moreover, by comparing the results of the two runs, one 
can see clearly that the cloud morphology is sensitive to the 
details of the thermal treatment adopted. In the run with 
the Koyama & Inutsuka cooling function, the clumps and 
filaments are contained within a region that is still roughly 
circular, and that has a radius that is only slightly larger 
than the initial radius of our inflowing gas (7? ~ 37 pc at 
t — 15 Myr, compared with R = 32 pc initially). On the 
other hand, in the non-equilibrium chemistry run, the dense 
gas occupies a significantly larger region, and is less circular, 
having an edge that is dominated by long, thin filaments of 
gas. The larger size of the cloud in the non-equilibrium run 
is also clearly apparent in the side-on view (Figure [SJ). 

Looking at the velocity distribution of the gas perpen- 
dicular to the flow (i.e. the radial velocities shown in Figure[6] 
for a slice through the centre of the dense cloud), we see that 
the immediate cause of the difference in morphologies is a 
difference in the velocity distributions. In the Koyama & In- 
utsuka run, the net outward velocity of the dense gas is very 
small. In the non-equilibrium chemistry run, on the other 
hand, the gas near the axis of the inflow is relatively static, 
but the dense gas close to the edges of the distribution is 
largely flowing outwards. In particular, the gas in the fila- 
ments has outward velocities of as much as 2 km s _1 . 

As the dense gas moves outwards in the non-equilibrium 
chemistry run, it drags the magnetic field lines along with it. 
The magnetic tension generated by this disturbance of the 
field lines exerts an inwards force on the expanding gas, and 
this force eventually becomes strong enough to reverse the 
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Figure 4. Column density of the inner region of the dense cloud viewed face-on at three different output times: t ~ 15 Myr, t ~ 25 Myr 
and t ~ 32 Myr. The top panels show the results from the run that used the Koyama & Inutsuka cooling function, while the bottom 
panels show the results from the run with the full non-equilibrium treatment. Note that output times considered here and in Figures [5HB1 
below are slightly different in the two runs owing to minor differences in the timing of the output snapshots produced by FLASH in the 
two different simulations. 




x [pc] x [pc] x [pc] 

Figure 5. Column density of the inner region of the dense cloud viewed edge-on at three different output times: t ~ 15 Myr, t ~ 25 Myr 
and t ~ 32 Myr. The top panels show the results from the run that used the Koyama & Inutsuka cooling function, while the bottom 
panels show the results from the run with the full non-equilibrium treatment. 
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Figure 6. Radial velocity of the gas in a slice through the centre of the cloud, relative to the gas at the centre of the cloud. As in 
Figures [4] and [5] we plot results from three different output times (i ~ 15 Myr, t ~ 25 Myr and t ~ 32 Myr) for both the simulation 
using the Koyama & Inutsuka cooling function (top panels) and the full non-equilibrium treatment (bottom panels). 



direction of the flow. This effect can be seen quite clearly 
when we look at the column density distribution and veloc- 
ity field at t ~ 25 Myr (the middle panels of Figures [4}[6]) . 
The gas distribution has become more compact, resulting 
in higher column densities, particularly close to the central 
axis of the flow, and more of the gas is flowing in than is 
flowing out. At an even later time (t ~ 32 Myr; right-hand 
panels in Figures [4j{6j) , the flow has "bounced" and has be- 
gun to re-expand once more. Looking at the results from the 
run with the Koyama & Inutsuka cooling function, we see 
hints of similar behaviour, but in this case both the initial 
outflow and the subsequent inflow are much weaker. 

The root cause of the difference in behaviour between 
the two runs is the thermal evolution of the shock-heated 
gas in the central cloud. The two inflowing streams of gas 
each are moving at a speed of 7 km s _1 , and so their rel- 
ative velocity is 14kms _1 , or around 2.5 times the speed 
of sound in the warm gas. When the gas collides, it passes 
through a shock, which heats it and compresses it. We can 
estimate the post-shock density and temperature by apply- 
ing the standard shock jump conditions. Since the magnetic 
field in our simulations is oriented along the flow, it plays no 
role in determining the post-shock conditions, and the same 
conditions apply as for a purely hydrodynamical shock. For 
the density, we therefore have the relationship 



p2 
pi 



(7- l)M 2 + 2' 



(11) 



ature, we have the relationship 
Ti 



(12) 



where p\ is the pre-shock density, p2 is the post-shock den- 
sity, and 7 is the adiabatic index of the gas. For the temper- 



where Ti and T 2 are the pre-shock and post-shock temper- 
atures, respectively. For the case of M = 2.5 and 7 = 5/3, 
we therefore find that p 2 ~ 2.7 pi and T 2 ~ 2.8 Ti. The 
post-shock thermal pressure is therefore a factor of around 
7.6 larger than the pre-shock thermal pressure, which itself 
is the same as the thermal pressure of the surrounding gas 
not participating in the inflow. The shocked gas is therefore 
over-pressured relative to its surroundings, and the resulting 
pressure gradient causes the gas to expand in the directions 
perpendicular to the inflow. In the inflow direction, it is 
confined by the ram pressure of the flow. The effectiveness 
with which this pressure gradient can accelerate gas out- 
wards from the central cloud depends on the length of time 
the gas remains in this over-pressured state. If we assume 
that the initial cooling of the gas is isochoric (i.e. that there 
is no change in its density), then it will cease to be over- 
pressured once its temperature drops below T\p\j p 2 , which 
for the case considered above and an initial temperature of 
5000 K yields T ~ 1850 K. 

In Figure [7] we show how the temperature of a par- 
cel of gas with density n — 2.7 cm -3 and temperature 
T = 14000 K changes as a function of time when we model 
heating and cooling using the Koyama & Inutsuka cooling 
function (dashed line) or our full non-equilibrium treatment 
(solid line). In the latter case, we take the initial chemical 
state of the gas to be the same as in our simulations. We 
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Figure 7. Evolution with time of the temperature of a fluid ele- 
ment with initial density rii = 2.7 cm -3 and initial temperature 
Tj = 14000 K, evolving at constant density. The initial chemical 
state of the gas, the strength of the UV radiation field, the cosmic 
ray ionisation rate and the metallicity are all taken to be the same 
as in our colliding flow simulations. The solid line shows the tem- 
perature evolution that we obtain when we model the gas using 
our non-equilibrium chemical model and cooling function, while 
the dashed line shows the results that we obtain when we use the 
Koyama & Inutsuka cooling function. The horizontal dash-dotted 
line shows the temperature at which the thermal pressure of the 
gas is the same as that of the unperturbed WNM in our colliding 
flow simulations. 

see that when we use the Koyama & Inutsuka cooling func- 
tion, the gas cools very rapidly at early times, returning to 
its original temperature after only 0.3 Myr. Subsequently, it 
cools more slowly, and it ceases to be over-pressured with 
respect to the unperturbed gas after around 1.3 Myr. On 
the other hand, when we model the chemistry and cooling 
of the gas using our non-equilibrium model, we see that it 
takes considerably longer for the gas too cool. In this case, 
the temperature of the gas returns to its original value after 
around 1.2 Myr, and the gas remains over-pressured until 
t ~ 2.45 Myr. These results demonstrate that when we use 
our non-equilibrium treatment to model the cooling of the 
gas, it remains significantly over-pressured for a much longer 
period than when we use the Koyama & Inutsuka cooling 
function. In the former case, therefore, the gas is acceler- 
ated by an outward-pointing pressure gradient for a longer 
period of time, and hence attains a significantly higher out- 
ward velocity, accounting for the differences we see in the 
morphology of the cloud and the velocity structure of the 
gas. 

3.3 Clump properties 

Finally, we explore whether the differences in cloud mor- 
phology and in the velocity field of the gas lead to significant 
differences in the statistical properties of the dense clumps 
formed in the flow. Although the fact that we do not include 
self-gravity in our models prevents us from following the fur- 
ther evolution of this dense gas in detail, we know from pre- 
vious studies that it is the ongoing growth and merger of 
these dense clumps that eventually leads to the formation 
of gravitation ally unstable pre-ste l lar cores and, ultimately, 
stars (see e.g. lBaneriee et al]|2009l ; IVazauez-Semadeni et al.l 



l201ll ). Differences in the clump properties at an early stage 
may therefore be indicative of differences in the ability of 
the clouds to form stars. 

We identify clumps in our simulations by searching for 
connected regions with densities above 50 cm -3 using the 
SEARCH3D routine of IDL. This routine is consecutively 
applied on the peak densities in the simulation excluding the 
already extracted clumps until we identified all connected 
regions with densities above 50 cm -3 . For reasons of com- 
putational efficiency, we restrict our search for clumps to 
radial distances R ^ 40 pc from the central axis of the flow, 
but from Figure [4] we can clearly see that this region con- 
tains almost the entire mass of dense gas, and so we are 
unlikely to miss many clumps. In Figure [5] we show some 
of the averaged internal properties of the clumps at times 
t ~ 15 Myr, t ~ 25 Myr, and t ~ 32 Myr for both models. 

The top panels of these figures show the ratio of 
the clump mass M to the local Jeans mass Mj, plot- 
ted as a function of the clump mass. The masses of the 
clumps span a wide range, 0.1-10 Mq. The fact that we 
find no clumps smaller than around 0.1 Mq is a conse- 
quence of the limited numerical resolution of our simula- 
tions, which suppresses the formation of extremely small 
scale structures. However, we see that even with our cur- 
rent resolution, the vast majority of the clumps are not 
self-gravitating, since they have M/Mj < 1. The only ex- 
ceptions are a couple of clumps present at t = 15 Myr 
that have M ~ 1-2 Mj. We therefore see that in order 
to form a significant number of self-gravitating clumps - 
a necessary pre-requisite for star formation - the compres- 
sions produced by the collision of the flows and the con- 
sequent thermal instability are not sufficient; some form of 
large-scale collapse of the cloud is also required. This is also 
evident in simulations including self-gravity, where in the 
case of a critical setu p the presence of magnetic field s sup- 
press star formation (|Vazquez-Semadeni et al1l201ll ). Sim- 
ilar results have been found in a number of other stud- 
ies (see e.g. iKovama fc Inutsukal 120021 : iHeitsch et alj 120051 ; 
IVazcmez-Semadeni et al.ll2007l ). 

Comparing the clump properties in the different mod- 
els, we see that the mean density of the clumps is not par- 
ticularly sensitive to the way in which the cooling of the 
gas is modelled, although there is a tentative hint that the 
clumps in the simulation with the Koyama & Inutsuka cool- 
ing function may be slightly denser than in those in the 
other simulation. A more pronounced difference is apparent 
in the mean temperature of the clumps, which we find is 
roughly 10 K higher when we use the non-equilibrium treat- 
ment than when we use the Koyama & Inutsuka cooling 
function. This difference is due to the difference in the pho- 
toelectric heating rate assumed in the two simulations, as 
noted in Section 13.11 above, and leads to minor differences 
in the Jeans masses that we determine for our clumps, as 
one can see from the upper panels in Figure [8] However, it 
should be noted that both simulations likely overestimate the 
temperature of these dense clump s, owing to their negle ct 
of the effects of d ust shielding (c.f. iGlover fc Clarldl2012l bl; 
IClark et alj|2012h . 




80 



60 



40 
1000 



100 



^ Koyama & Inutsuka 
O Glover & Mac Low 



I II 




I I I M ill — I I I M ill 



-+++++1 — I II 




10 

10' 

S 10" 3 

lO" 4 

lO" 5 

lO" 6 
100 



80 

'0 



60- 



40 
1000 



0.1 



10 100 1000 10000 
M [M ] 



100 r 



0.1 



25 Myr 





8* <s> °* * * 




10 100 1000 10000 
M [M ] 




10 100 1000 10000 
M [M ] 



Figure 8. Properties of the set of dense clumps identified as described in Section [3.3l Results are shown for three output times: t ~ 15 Myr 
(left-hand panels), t ~ 25 Myr (central panels) and t ~ 32 Myr (right-hand panels). The upper row of panels shows the ratio of the clump 
mass to the local Jeans mass, Mj, the central row shows the mean density of the clumps and the lower row shows the mean temperature. 
In each case, we plot these quantities as a function of the total clump mass. The red symbols indicate the results that we obtain in the 
run with our non-equilibrium chemical model, while the black symbols indicate the results that we obtain when we use the Koyama & 
Inutsuka cooling function. 



4 SUMMARY 

We have presented the results of a study that examines the 
influence of two different thermal models on the formation 
of cold, dense clouds within converging flows of warm atomic 
gas, and on the nature of the clumps that form within these 
clouds. To do this, we performed high-resolution 3D MHD 
simulations using the massively parallel code FLASH, mod- 
ified to include a detailed treatment of atomic and molecular 
cooling, and a simplified but acc urate treatment of the most 
important hydrogen chemistry (|Micic et al.ll2012T ). We di- 
rectly compare the results obtained from this model with 
those that we obtain i f we use a simplified cooling function, 
taken from the work of lKovama fc Inutsukal (2002), that has 
been used in a number of o ther studies of cloud formation 
in converging flows (see e.g. iBaneriee et al]|2009f ). 

We find that the density and temperature PDFs pro- 
duced in the two simulations are qualitatively similar, al- 
though some minor quantit ative differences exist. In com- 
mon with previous work (e.g.lVazauez-Semadeni et al ]|2007l : 



iHennebelle et al.ll200sl ; IBaneriee et al.ll2009h we find that in 
addition to clear CNM and WNM phases, there is also a sig- 



nificant fraction of mass located in the thermally unstable 
region between these two phases. The temperatures recov- 
ered for the CNM and WNM phases are slightly smaller 
in the simulation run with the Koyama & Inutsuka cooling 
function, as their approach significantly overestimates the 
cooling rate of hot (T ~ 10 4 K) gas and underestimates the 
photoelectric heating rate in the CNM regime by a factor 
of 2.5. In the Koyama & Inutsuka model, most of the gas 
is in thermal equilibrium, while in the full non-equilibrium 
model, significant deviations from thermal equilibrium are 
apparent for gas in the density range 1 < n < 30 cm -3 . In 
this intermediate density regime, the photoelectric heating 
efficiency has not yet reached its limiting value and increases 
with increasing electron number density. Due to the rela- 
tively long recombination timescale, the gas is slightly more 
ionised than it would be in chemical equilibrium, and hence 
is heated slightly more efficiently. 

We have also shown that the cloud morphology is sen- 
sitive to the choice of the thermal treatment. The cloud 
formed in the non-equilibrium chemistry run is larger and 
more filamentary than the cloud formed in the run with the 
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Koyama & Inutsuka cooling function. This difference in mor- 
phology is caused by a difference in the velocity distribution 
of the gas, which itself can be understood as a consequence 
of the difference in the cooling time of hot gas in the two 
models. In the Koyama & Inutsuka model, the high temper- 
ature, shock-heated gas cools very rapidly, and the gas re- 
mains over-pressure with respect to its surroundings for only 
a very short time. With the full non-equilibrium treatment, 
however, the cooling time is significantly longer and hence 
the cloud remains over-pressured for longer. Consequently, 
the gas is accelerated outwards more efficiently, resulting in 
a larger, more disordered gas distribution. 

Finally, we have investigated whether the properties of 
the clumps that form in the clouds differ significantly be- 
tween the two runs. We have identified clumps with masses 
in a range of 0.1-10 3 Mq in both simulations, but find that 
almost all of these structures are not self-gravitating, having 
M/Mj < 1. This suggests that the compressions produced 
by the collision of the flows followed by thermal instability 
are not sufficient to lead to the formation of gravitationally 
unstable pre-stellar cores and stars. In order for star for- 
mation to take place, some form of large-scale collapse or 
compression of the cloud in a direction perpendicular to the 
inflow appears to be required. One possible source for this 
compression is the global gravitation al contraction of the 
shocked cloud, as previously note d in [K oyama & Inutsuka 
1 20021) . iHeitsch et all (|2005l ). and IVazq uez-Semadeni et aJj 
12007). Another possibility is compression due to the ef- 
fects of orbit crossing in a spira l arm shock, as discussed in 
iBonnell. Dobbs fc Smith! (|2013T ). Most of the properties of 
the clumps are very similar in our two models, with the most 
significant difference being a systematic offset of roughly 
10 K between the mean clump temperatures in the Koyama 
& Inutsuka model and those in the non-equilibrium model. 
However, it should be noted that both models overestimate 
the temperature of the clumps, owing to their neglect of the 
effects of dust shielding. 
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